#include "ms_voxelimage.h"
#include "fe_vecmat.h"
#include "ms_misc.h"
#include <sys/stat.h>
#include <iostream>
#include <fstream>

// #include <ImageMagick/magick/MagickCore.h>
#include <magick/MagickCore.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#ifndef PI
#define UNDEF_PI
#define PI 3.1415926535897932384626433
#endif



using namespace std;


// class VoxelImage {
// 	public:
// 		unsigned char *** data;
// 		int Nx, Ny, Nz;
// 		const double h;
// 
// 	public:
// 		VoxelImage( const double L = 1., const int N = 10 );
// 		VoxelImage( const double Lx, const double Ly, const double Lz, const double h );
// 
// 		void	Erase( const unsigned char bgrd_col = 255 );
// 		void	DrawBox( const double x0[3], const double x1[3], const unsigned char col = 0 );
// 		void	DrawCircle( const double x[0], const double r, const unsigned char col = 0 );
// 		void	WriteImageStack( const char fn[] );
// 		void	WriteVTKFile( const char fn[] );
// };

/* ******************************************************************* */
VoxelImage::VoxelImage( const double L, const int N, const unsigned char col )
{
  
	length = L;
/*	if(N % 2 != 1) n++;
	if(n<3) n = 3;
// 	h = L/(double)(n-1);*/
	wx = wy = wz = 1.;
	h = L/(double)N;
	Nx = Ny = Nz = N;

	data = new unsigned char ** [Nx];
	for(int ix = 0; ix < Nx; ix++)
	{
		data[ix] = new unsigned char * [Ny];
		for(int iy = 0; iy < Ny; iy++)
		{
			data[ix][iy] = new unsigned char [Nz];
			memset(data[ix][iy], col, Nz);
		}
	}
}
/* ******************************************************************* */
void VoxelImage::SetData(int x, int y, int z){

data[x][y][z] = true;

}
/* ******************************************************************* */
void VoxelImage::GetVoxelCoordinates( const double x, const double y, const double z, int & ix, int & iy, int & iz ) const
{
	ix = GetVoxelCoX(x);
	iy = GetVoxelCoX(y);
	iz = GetVoxelCoX(z);
}
/* ******************************************************************* */
VoxelImage::VoxelImage( const double Lx, const double Ly, const double Lz, const double H, const unsigned char col )
{
	h	= H;
	Nx	= (int)(Lx / h);	Ny	= (int)(Ly / h);	Nz	= (int)(Lz / h);
	wx	= H*Nx;				wy	= H*Ny;				wz	= H*Nz;

	data = new unsigned char ** [Nx];
	for(int ix = 0; ix < Nx; ix++)
	{
		data[ix] = new unsigned char * [Ny];
		for(int iy = 0; iy < Ny; iy++)
		{
			data[ix][iy] = new unsigned char [Nz];
			memset(data[ix][iy], col, Nz);
		}
	}
}
/* ******************************************************************* */
VoxelImage::~VoxelImage()
{
	for(int ix=0; ix<Nx; ix++) {
		for(int iy=0; iy<Ny; iy++) {
			delete [] data[ix][iy]; data[ix][iy] = 0;
		}
		delete [] data [ix]; data[ix] = 0;
	}
	delete [] data; data = 0;
}
/* ******************************************************************* */
void	VoxelImage::Erase( const unsigned char bgrd_col  ) {
	for(int ix = 0; ix < Nx; ix++)
		for(int iy = 0; iy < Ny; iy++)
			memset(data[ix][iy], bgrd_col, Nz);
}
/* ******************************************************************* */
void	VoxelImage::DrawBox( const double x0[3], const double x1[3], const unsigned char col ) {

	int ix0,iy0,iz0, ix1,iy1,iz1;

	GetVoxelCoordinates( x0[0],x0[1],x0[2], ix0, iy0, iz0 );
	GetVoxelCoordinates( x1[0],x1[1],x1[2], ix1, iy1, iz1 );

	if( ix0 > ix1 ) { int i = ix1; ix1=ix0; ix0=i; }
	if( iy0 > ix1 ) { int i = iy1; iy1=iy0; iy0=i; }
	if( iz0 > ix1 ) { int i = iz1; iz1=iz0; iz0=i; }
	if( ix0 >= Nx) return;
	if( ix1 < 0 ) return;
	if( iy0 >= Ny) return;
	if( iy1 < 0 ) return;
	if( iz0 >= Ny) return;
	if( iz1 < 0 ) return;

	if( ix0 < 0 ) ix0  = 0;
	if( ix1 >= Nx) ix1 = Nx-1;
	if( iy0 < 0 ) iy0  = 0;
	if( iy1 >= Ny) iy1 = Ny-1;
	if( iz0 < 0 ) iz0  = 0;
	if( iz1 >= Nz) iz1 = Nz-1;
	for(int ix=ix0; ix<=ix1; ix++)
		for(int iy=iy0; iy<=iy1; iy++)
			memset( data[ix][iy] + iz0, col, iz1 - iz0 + 1 );
}
/* ******************************************************************* */
int		VoxelImage::DrawCircleZ( const double x[], const double r, const unsigned char col ) {
	int		ix0, iy0, iz0;
	GetVoxelCoordinates( x[0], x[1], x[2], ix0, iy0, iz0 );

	if( ( iz0 <  0  ) || ( iz0 >= Nz ) ) {
		fprintf(stderr,"ERROR: iz0 = %5d   (  x  =  [  %16.8lf,  %16.8lf,  %16.8lf  ]  )\n", iz0, x[0], x[1], x[2] );
		return 1;
	}

	double	x0	= GetX(ix0),
			y0	= GetY(iy0);
	
	int ix_start	= GetVoxelCoX( x[0]-r ),
		ix_end		= GetVoxelCoX( x[0]+r ),
		iy_start, iy_end;

	if( ( ix_end <  0  ) || ( ix_start >= Nx ) ) {
/*		fprintf(stderr,"ERROR: ix_start = %5d ( Nx = %5d )   ix_end = %5d ( Nx = %5d )   (  x  =  [  %16.8lf,  %16.8lf,  %16.8lf  ]  )\n",
				ix_start, Nx, ix_end, Nx, x[0], x[1], x[2] );
		fprintf(stderr,"r=%16.8lf, (int)((x[0]-r)/h) = %d  + int((Nx-1)/2) =%d \n", r, (int)((x[0]-r)/h) , int((Nx-1)/2));
		fprintf(stderr,"h=%16.8lf\n", h); */
		return 1;
	}
	
	const int	iy_0	= GetVoxelCoY( x[1]-r ),
				iy_1	= GetVoxelCoY( x[1]+r );

	if( ( iy_0 >= Ny ) || ( iy_1 < 0 ) ) {
		// fprintf(stderr,"ERROR: iy0 = %5d    iy1 = %5d   (  x  =  [  %16.8lf,  %16.8lf,  %16.8lf  ]  )\n", iy_0, iy_1, x[0], x[1], x[2] );
		return 1;
	}

	if( ix_start <  0 ) ix_start = 0;
	if( ix_end   >= Nx) ix_end   = Nx-1;
	
	double X,by;
	int n;
	
	X = GetX(ix_start);
	for(int ix=ix_start; ix<=ix_end; ix++)
	{
		by			= r*r-(X-x[0])*(X-x[0]);
		if(by < 0.)		by	= 0.;
		else			by	= sqrt(by);
		iy_start	= GetVoxelCoY(x[1]-by);
		iy_end		= GetVoxelCoY(x[1]+by);
		
		if( (iy_start < Ny) && (iy_end >= 0) )
		{
			if( iy_start <  0  ) iy_start = 0;
			if( iy_end   >= Ny ) iy_end   = Ny-1;
			for( int iy=iy_start; iy<=iy_end; iy++ )
				data[ix][iy][iz0] = col;
		}
		X += h;
	}
	return 0;
}
/* ******************************************************************* */
int VoxelImage::ShowNumber( char fn[] )
/** 
  * Diese Funktion erzeugt die charakteristischen Zahlen des Spongebobs
  * @Param kein Parameter wird übergeben
  * @Return Alle charakteristischen Zahlen werden normal mit cout ausgegeben, 
	    danach wird die Gesamtsumme der charakteristischen Zahl übergeben.
 */

{
	// Parameter initialisieren 
	int sum  = 0; // Summe des Histogramms
	int sum2 = 0; // Gesamtsumme der binären charakteristischen Zahl
	int sum3 = 0; // Summe aller Histogramm-Werte
	int sum4 = 0; // Summe aller Histogramm-Werte, in denen der erste Punkt belegt ist (ungerade)
	
	double volfrac = 0.; // volume fraction
	double specsurf = 0.; // specific surface area
	double meancurv = 0.; // Mean curvature
	double totalcurv = 0.; // Total curvature
	double eulerzahl = 0.; // Euler Zahl
	double a[13];
	
	int delta[3]; // Kantenlaenge des Wuerfels
	int hist[256]; // Anzahl der Varianten der char. Zahl
	
	// Parameter fuer spezifische Oberflaeche
	int koeff[13][2] = {
          {1,2},{1,4},{1,16}, //
          {1,8},{2,4},{1,32},{2,16},{1,64},{4,16}, // Diagonal Ebene
          {1,128},{2,64},{4,32},{8,16} // Diagonal Raum
	};
	
	// Parameter fuer mittlere Kruemmung
        int krue[9][4] = {
        {1,2,4,8} , {1,2,16,32} , {1,4,16,64} , {1,2,64,128} , {4,16,8,32} ,
        {1,32,4,128} , {2,8,16,64} , {2,4,32,64} , {1,16,8,128}
	};
        int kte[8][3] = {
          {1,64,32} , {2,16,128} , {8,64,32} , {4,16,128} ,
          {2,4,128} , {8,1,128} , {2,4,16} , {8,1,32}
	};
	
	
	double cv[13] = { // Abstaende der n'ten Voronoi-Zelle, fuer: Mittlere Kruemmung, SpecSurfaceArea
	0.045778 , 0.045778 , 0.045778 ,
	0.036981 , 0.036981 , 0.036981 , 0.036981 , 0.036981 , 0.036981 ,
	0.035196 , 0.035196 , 0.035196 , 0.035196
	};
	
	// befuellen Delta mit 3 Werten = Kantenlaenge
	for (int i=0;i<3;i++){
		delta[i]=length;
	}
	
	double r[13]; // Abstaende
	r[0] = r[1] = r[2] = delta[0]; // Kantenlaenge
	r[3] = r[4] = r[5] = r[6] = r[7] = r[8] = sqrt(delta[0]*delta[0] + delta[1]*delta[1] ); // Abstaende Diagonalpunkte Ebene
	r[9] = r[10] =r[11] = r[12] = sqrt(delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2] ); // Abstaende Diagonalpkte Volumen 
	
	// fuer Mittlere Kruemmung
	double si = (r[0]+r[1]+r[3])/2; // Heron Formula
	
	// 
	a[0] = delta[0]*delta[1];
	a[1] = delta[0]*delta[2];
	a[2] = delta[1]*delta[2];
	a[3]= a[4] = a[5] = a[6] = a[7] = a[8] = delta[2]*r[3];
	a[9]= a[10] = a[11] = a[12] = 2 * sqrt(si*(si-r[3])*(si-r[4])*(si-r[5])); // Triangelflaeche


	
	// befuellen des Histogramms mit 256 Werten
	for (int i=0;i<256;i++){
		hist[i]=0;
	}

	
	// erstellen der charakteristischen Zahl
	for(int ix=0;ix<Nx;ix++)
	{
		for(int iy=0;iy<Ny;iy++) 
		{
			for(int iz=0;iz<Nz;iz++)
			{
				if(data[ix]		[iy]		[iz]) 		sum += 1; // wenn dat = true -> summe erhoehen
				if(data[(ix+1) % Nx]	[iy]		[iz]) 		sum += 2;
//				if(data[(ix+1) % Nx]	[(iy+1) % Ny]	[iz]) 		sum += 4;
                                if(data[ix]		[(iy+1) % Ny]	[iz]) 		sum += 4;
                                if(data[(ix+1) % Nx]	[(iy+1) % Ny]	[iz]) 		sum += 8;


				if(data[ix]		[iy]		[(iz+1) % Nz]) 	sum += 16;
				if(data[(ix+1) % Nx]	[iy]		[(iz+1) % Nz]) 	sum += 32;
//				if(data[(ix+1) % Nx]	[(iy+1) % Ny]	[(iz+1) % Nz]) 	sum += 64;
                                if(data[ix]		[(iy+1) % Ny]	[(iz+1) % Nz]) 	sum += 64;
                                if(data[(ix+1) % Nx]	[(iy+1) % Ny]	[(iz+1) % Nz]) 	sum += 128;


				sum2 = sum2 + sum; //Gesamtsumme erhöhen
				hist[sum] = hist[sum]+1; // Histogramm an der richtigen Stelle eines erhoehen
				
				
				sum = 0; // Summe wieder loeschen fuer naechsten Summe
				
			}
		}
	}

	// Volumendichte, lokale Volumenfraktion:
        for(int l=0;l<256;l++) sum3 = sum3 + hist[l];
	for(int k=0;k<128;k++) sum4 = sum4 + hist[2*k + 1];
	volfrac = static_cast<float>(sum4)/static_cast<float>(sum3);

	// spezifische Oberfläche

        specsurf = 0.;
	for(int i=0;i<256;i++){
		for(int a=0;a<13;a++){
			specsurf += hist[i] * cv[a]/r[a] * ((i==(i|koeff[a][0]))*(0==(i&koeff[a][1])) + (i==(i|koeff[a][1]))*(0==(i&koeff[a][0])));
		}

	}

        specsurf = 4.0*specsurf/static_cast<float>(sum3);
	
	// Mittlere Kruemmung
	for(int l=0; l<256;l++)
	{
	  for(int kx=0;kx<9;kx++){
	    for(int ky=0;ky<4;ky++)
	    {
	      meancurv += static_cast<double>(hist[l])*cv[kx] / (4.0*a[kx]) 
	      * ((l==(l|krue[kx][ky]))*(0==(l&krue[kx][(ky+1)%4]))
	      * (0==(l&krue[kx][(ky+2)%4]))*(0==(l&krue[kx][(ky+3)%4]))
	      - (l==(l|krue[kx][(ky)]))*(l==(l|krue[kx][(ky+1)%4]))
	      * (l==(l|krue[kx][(ky+2)%4]))*(0==(l&krue[kx][(ky+3)%4])));
	    } 
	  }    
	  for(int kx= 9;kx<13;kx++)
	  {
	    for(int ky=0;ky<3;ky++)
	    {
	      meancurv += static_cast<double>(hist[l])*cv[kx] / (4.0*a[kx]) 
	      * ((l==(l|kte[kx-9][ky]))*(0==(l&kte[kx-9][(ky+1)%3]))
	      * (0==(l&kte[kx-9][(ky+2)%3]))
	      - (l==(l|kte[kx-5][(ky)]))*(l==(l|kte[kx-5][(ky+1)%3]))
	      * (0==(l&kte[kx-5][(ky+2)%3])));
	    } 
	  }
	}
	meancurv = 4.0 * PI * meancurv / sum3;
	
	// Euler-Zahl
	
	long int Eul = 0;
	int euler[256] = {
	  0 , 1 , 0 , 0 , 0 , 0 , 0 ,-1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,
	  0 , 0 , 0 ,-1 , 0 ,-1 , 0 ,-2 , 0 , 0 , 0 ,-1 , 0 ,-1 , 0 ,-1 ,
	  0 , 1 , 0 , 0 , 0 , 0 , 0 ,-1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,
	  0 , 0 , 0 , 0 , 0 ,-1 , 0 ,-1 , 0 , 0 , 0 , 0 , 0 ,-1 , 0 , 0 ,
	  0 , 1 , 0 , 0 , 0 , 0 , 0 ,-1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,
	  0 , 0 , 0 ,-1 , 0 , 0 , 0 ,-1 , 0 , 0 , 0 ,-1 , 0 , 0 , 0 , 0 ,
	  0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,
	  0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 ,
	  0 , 1 , 0 , 0 , 0 , 0 , 0 ,-1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,
	  0 , 0 , 0 ,-1 , 0 ,-1 , 0 ,-2 , 0 , 0 , 0 ,-1 , 0 ,-1 , 0 ,-1 ,
	  0 , 1 , 0 , 0 , 0 , 0 , 0 ,-1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,
	  0 , 0 , 0 , 0 , 0 ,-1 , 0 ,-1 , 0 , 0 , 0 , 0 , 0 ,-1 , 0 , 0 ,
	  0 , 1 , 0 , 0 , 0 , 0 , 0 ,-1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,
	  0 , 0 , 0 ,-1 , 0 , 0 , 0 ,-1 , 0 , 0 , 0 ,-1 , 0 , 0 , 0 , 0 ,
	  0 , 1 , 0 , 0 , 0 , 0 , 0,  0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 ,
	  0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0
	};
	for(int l=0; l<256;l++) Eul += euler[l]*hist[l];
	eulerzahl = static_cast<double>(Eul) / static_cast<double>(sum3);


	
	// Totale Kruemmung
	
	long int Chi = 0;
	int iv[256] = {
	  0 , 3 , 3 , 0 , 3 , 0 , 6 ,-3 , 3 , 6 , 0 ,-3 , 0 ,-3 ,-3 , 0 ,
	  3 , 0 , 6 ,-3 , 6 ,-3 , 9 ,-6 , 6 , 3 , 3 ,-6 , 3 ,-6 , 0 ,-3 ,
	  3 , 6 , 0 ,-3 , 6 , 3 , 3 ,-6 , 6 , 9 ,-3 ,-6 , 3 , 0 ,-6 ,-3 ,
	  0 ,-3 ,-3 , 0 , 3 ,-6 , 0 ,-3 , 3 , 0 ,-6 ,-3 , 0 ,-8 ,-8 , 0 ,
	  3 , 6 , 6 , 3 , 0 ,-3 , 3 ,-6 , 6 , 9 , 3 , 0 ,-3 ,-6 ,-6 ,-3 ,
	  0 ,-3 , 3 ,-6 ,-3 , 0 , 0 ,-3 , 3 , 0 , 0 ,-8 ,-6 ,-3 ,-8 , 0 ,
	  6 , 9 , 3 , 0 , 3 , 0 , 0 ,-8 , 9 ,12 , 0 ,-3 , 0 ,-3 ,-8 ,-6 ,
	 -3 ,-6 ,-6 ,-3 ,-6 ,-3 ,-8 , 0 , 0 ,-3 ,-8 ,-6 ,-8 ,-6 ,-12, 3 ,
	  3 , 6 , 6 , 3 , 6 , 3 , 9 , 0 , 0 , 3 ,-3 ,-6 ,-3 ,-6 ,-6 ,-3 ,
	  6 , 3 , 9 , 0 , 9 , 0 ,12 ,-3 , 3 , 0 , 0 ,-8 , 0 ,-8 ,-3 ,-6 ,
	  0 , 3 ,-3 ,-6 , 3 , 0 , 0 ,-8 ,-3 , 0 , 0 ,-3 ,-6 ,-8 ,-3 , 0 ,
	 -3 ,-6 ,-6 ,-3 , 0 ,-8 ,-3 ,-6 ,-6 ,-8 ,-3 , 0 ,-8,-12 ,-6 , 3 ,
	  0 , 3 , 3 , 0 ,-3 ,-6 , 0 ,-8 ,-3 , 0 ,-6 ,-8 , 0 ,-3 ,-3 , 0 ,
	 -3 ,-6 , 0 ,-8 ,-6 ,-3 ,-3 ,-6 ,-6 ,-8 ,-8,-12 ,-3 , 0 ,-6 , 3 ,
	 -3 , 0 ,-6 ,-8 ,-6 ,-8 ,-8,-12 ,-6 ,-3 ,-3 ,-6 ,-3 ,-6 , 0 , 3 ,
	  0 ,-3 ,-3 , 0 ,-3 , 0 ,-6 , 3 ,-3 ,-6 , 0 , 3 , 0 , 3 , 3 , 0
	};
	for(int l=0; l<256;l++) Chi += iv[l]*hist[l];
	totalcurv = PI/6. * static_cast<double>(Chi) / static_cast<double>(sum3)*delta[0]*delta[1]*delta[2];
	
  
	
	
	// strcat(fn , "/Data.csv");

	// Einfügen der Daten in eine .csv Datei
	ofstream scriptfile;
	scriptfile.open( fn , ios::app|ios::out );
        scriptfile << sum2 << ";" << volfrac << ";" << specsurf << ";" << meancurv << ";" << eulerzahl << ";" << totalcurv << endl;
	
	scriptfile.close();

	
	
	// Werte ausspucken:
	// for(int i=0;i<256;i++) cout << "Histogramm: " << i << " - "<< hist[i] << endl;
	cout << "charakteristische Zahl: " << sum2 << endl;
	cout << "Volumenfraktion: " << volfrac << endl;
	cout << "Spezifische Volumenfläche: " << specsurf << endl;
	cout << "Mittlere Krümmung: " << meancurv << endl;
	cout << "Eulerzahl: " << eulerzahl << endl;
	cout << "Gesamtkrümmung: " << totalcurv << endl;

	// cout << "sum3: " << sum3 << endl;
	// cout << "sum4: " << sum4 << endl;

	return sum2;
}



/* ******************************************************************* */
void	VoxelImage::WriteImageStack( const char fn[] ) const {

	/* new image pointer */
	Image			*img	= 0;
	
	/* single slice of 3d image: 3 rgb values for each pixel */
	float			*px		= (float*)malloc( 3*sizeof(float)*Nx*Ny );
	
	/* catch error messages */
	ExceptionInfo	*exc;
	
	/* some variables */
	char			fullname[1024];
	int ct=0;
	
	/* ************************************************ */
	MagickCoreGenesis(".", MagickTrue );
	ImageInfo	*info = CloneImageInfo((ImageInfo *) NULL);

	exc		= AcquireExceptionInfo();

	// printf("* Image resolution (nx, ny, nz): %d, %d, %d \n",Nx,Ny,Nz);

	for(int iz=0;iz<Nz;iz++)
	{
	  ct=0;
	 // printf("* Process slice %d of %d\n",iz,Nz);
		for(int iy=0;iy<Ny;iy++) 
			for(int ix=0;ix<Nx;ix++)
			{
				px[ct] = px[ct+1] = px[ct+2] = (data[ix][iy][iz]!=0)?1.:0.;
				ct+=3;
			}
	
		/* build one image and write to PNG file */
		img = ConstituteImage(
				Nx, Ny,  "RGB", FloatPixel,
				px, exc
			);
			
		sprintf(fullname,"%s%03d.png", fn,iz+1);
		strcpy( img->filename, fullname );
// 		img->compression = NoCompression;
		WriteImage( info, img );
		img=DestroyImage(img);

//		z += dx;
	}
	MagickCoreTerminus();
	free(px);
	char cmd[512];
	for(int i=0;i<Nz;i++) {
		sprintf(cmd,"convert -quality 100 %s%03d.png -resize 128x128 %s%03d.jpg && rm %s%03d.png", fn,i+1, fn,i+1, fn,i+1);
		system(cmd);
	}

}
/* ******************************************************************* */
void	VoxelImage::WriteVTKFile( const char fn[] ) const {
	FILE * F = fopen(fn,"wb");
	fprintf(F,"# vtk DataFile Version 2.0\n");
	fprintf(F,"Volume example\n");
// 	fprintf(F,"ASCII\n");
	fprintf(F,"BINARY\n");
	fprintf(F,"DATASET STRUCTURED_POINTS\n");
	fprintf(F,"DIMENSIONS %d %d %d\n", Nx, Ny, Nz);
	fprintf(F,"ASPECT_RATIO 1 1 1\n");
	fprintf(F,"ORIGIN %d %d %d\n", (Nx+1)/2, (Ny+1)/2, (Nz+1)/2);
	fprintf(F,"POINT_DATA %d\n", Nx*Ny*Nz);
	fprintf(F,"SCALARS material_number char 1\n");
	fprintf(F,"LOOKUP_TABLE default\n");
	//+++ unsigned char* xrow = new unsigned char[Nx]; ?
	unsigned char xrow[Nx];
	for(int iz=0;iz<Nz;iz++)
		for(int iy=0;iy<Ny;iy++) {
			for(int ix=0;ix<Nx;ix++)
				xrow[ix]=data[ix][iy][iz];
			fwrite(xrow, 1, Nx, F);
		}
	//+++ delete[] xrow; ?

	fprintf(F,"\n");
	fclose(F); F=0;
}

/* ******************************************************************* */
void VoxelImage::DrawLine(const int X0[], const int X1[], const unsigned short col )
{
	int		x0 = X0[0], x1 = X1[0],
			y0 = X0[1], y1 = X1[1],
			z0 = X0[2], z1 = X1[2];
	//'steep' xy Line, make longest delta x plane  
	
	const bool swap_xy = ( abs(y1-y0) > abs(x1-x0) );
	if ( swap_xy )
	{	Swap(x0, y0);	Swap(x1, y1); }

	const bool swap_xz = ( abs(z1-z0) > abs(x1-x0) );
	if ( swap_xz ) //+++
	{	Swap(x0, z0);	Swap(x1, z1); }
	
	int		dx = x1-x0,
			dy = y1-y0,
			dz = z1-z0;
	const int delta_x = abs(dx), delta_y = abs(dy), delta_z = abs(dz);
	
	//drift controls when to step in 'shallow' planes
	//starting value keeps Line centred
	int drift_xy = (delta_x / 2),
		drift_xz = (delta_x / 2),
		step_x   = 1,
		step_y   = 1,
		step_z   = 1;
	if (x0 > x1) step_x = -1;
	if (y0 > y1) step_y = -1;
	if (z0 > z1) step_z = -1;
	
	//starting point
	int x, y = y0, z = z0;
	int cx, cy, cz;
	for( x = x0; x != x1; x += step_x )
	{
		//copy position
		cx = x;		cy = y;		cz = z;

		//unswap (in reverse)
		if (swap_xz) Swap(cx, cz);
		if (swap_xy) Swap(cx, cy);
		
		data[cx][cy][cz] = col;
		
		//update progress in other planes
		drift_xy -= delta_y;
		drift_xz -= delta_z;

		//step in y plane
		if( drift_xy < 0 )
		{
			y			+= step_y;
			drift_xy	+= delta_x;
		}
		
		//same in z
		if( drift_xz < 0 )
		{
			z			+= step_z;
			drift_xz	+= delta_x;
		}
	}
	
}

/* ******************************************************************* */
void VoxelImage::DrawLine(const double X0[], const double X1[], const unsigned short col )
{
	double	x = X0[0],		y = X0[1],		z = X0[2],		L = Distance(X0, X1);
// 	const int N = L/h;
	const double	lx = fabs(X1[0]-x), ly = fabs(X1[1]-y), lz = fabs(X1[2]-z);
	const double	mx = lx/h, my = ly/h, mz = lz/h; //+++ double /int war hier ursprünglich
	
	const int N = int(  1.1 *( ( mx < my ) ? ( ( my < mz ) ? mz : my ) : ( ( mx < mz ) ? mz : mx ) ) ); 
	//+++ Die Werte mx, my, mz als int haben keine Nachkommastellen? 
	const double	dx = (X1[0]-x) / (double)N,
					dy = (X1[1]-y) / (double)N,
					dz = (X1[2]-z) / (double)N;
	int ix,iy,iz;
	bool sc = false;
// 	int ix_start,iy_start,iz_start;
// 	int ix_end,iy_end,iz_end;
// 	int i=0;
// 	for ( ; (i<=N) && !sc ; i++) {
// 		GetVoxelCoordinates( x, y, z, ix_start, iy_start, iz_start );
// 		sc = ( XAdmissible(ix_start) && YAdmissible(iy_start) && ZAdmissible(iz_start) );
// 		x += dx; y += dy, z += dz;
// 	}
// 	sc=false;
// 	for ( ; (i<=N) && !sc ; i++) {
// 		GetVoxelCoordinates( x, y, z, ix_end, iy_end, iz_end );
// 		sc = ( XAdmissible(ix_end) && YAdmissible(iy_end) && ZAdmissible(iz_end) );
// 		x += dx; y += dy, z += dz;
// 	}
// 	DrawLine(	(int []){ix_start, iy_start, iz_start }, 
// 				(int []){ix_end, iy_end, iz_end },col );
// 	
	
	for (int i=0; (i<=N) ; i++)
	{
		GetVoxelCoordinates( x, y, z, ix, iy, iz );
		sc = ( XAdmissible(ix) && YAdmissible(iy) && ZAdmissible(iz) );
		if(sc) {
			data[ix][iy][iz] = col;
// 			printf("%6.3lf, %6.3lf, %6.3lf\n",x,y,z);
		}
		x += dx; y += dy; z += dz;
	}
}
